<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>claqr5.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="CLAQR5.1"></a><a href="claqr5.f.html#CLAQR5.1">CLAQR5</a>( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S,
     $                   H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV,
     $                   WV, LDWV, NH, WH, LDWH )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK auxiliary routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      INTEGER            IHIZ, ILOZ, KACC22, KBOT, KTOP, LDH, LDU, LDV,
     $                   LDWH, LDWV, LDZ, N, NH, NSHFTS, NV
      LOGICAL            WANTT, WANTZ
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      COMPLEX            H( LDH, * ), S( * ), U( LDU, * ), V( LDV, * ),
     $                   WH( LDWH, * ), WV( LDWV, * ), Z( LDZ, * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     This auxiliary subroutine called by <a name="CLAQR0.19"></a><a href="claqr0.f.html#CLAQR0.1">CLAQR0</a> performs a
</span><span class="comment">*</span><span class="comment">     single small-bulge multi-shift QR sweep.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      WANTT  (input) logical scalar
</span><span class="comment">*</span><span class="comment">             WANTT = .true. if the triangular Schur factor
</span><span class="comment">*</span><span class="comment">             is being computed.  WANTT is set to .false. otherwise.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      WANTZ  (input) logical scalar
</span><span class="comment">*</span><span class="comment">             WANTZ = .true. if the unitary Schur factor is being
</span><span class="comment">*</span><span class="comment">             computed.  WANTZ is set to .false. otherwise.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      KACC22 (input) integer with value 0, 1, or 2.
</span><span class="comment">*</span><span class="comment">             Specifies the computation mode of far-from-diagonal
</span><span class="comment">*</span><span class="comment">             orthogonal updates.
</span><span class="comment">*</span><span class="comment">        = 0: <a name="CLAQR5.33"></a><a href="claqr5.f.html#CLAQR5.1">CLAQR5</a> does not accumulate reflections and does not
</span><span class="comment">*</span><span class="comment">             use matrix-matrix multiply to update far-from-diagonal
</span><span class="comment">*</span><span class="comment">             matrix entries.
</span><span class="comment">*</span><span class="comment">        = 1: <a name="CLAQR5.36"></a><a href="claqr5.f.html#CLAQR5.1">CLAQR5</a> accumulates reflections and uses matrix-matrix
</span><span class="comment">*</span><span class="comment">             multiply to update the far-from-diagonal matrix entries.
</span><span class="comment">*</span><span class="comment">        = 2: <a name="CLAQR5.38"></a><a href="claqr5.f.html#CLAQR5.1">CLAQR5</a> accumulates reflections, uses matrix-matrix
</span><span class="comment">*</span><span class="comment">             multiply to update the far-from-diagonal matrix entries,
</span><span class="comment">*</span><span class="comment">             and takes advantage of 2-by-2 block structure during
</span><span class="comment">*</span><span class="comment">             matrix multiplies.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      N      (input) integer scalar
</span><span class="comment">*</span><span class="comment">             N is the order of the Hessenberg matrix H upon which this
</span><span class="comment">*</span><span class="comment">             subroutine operates.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      KTOP   (input) integer scalar
</span><span class="comment">*</span><span class="comment">      KBOT   (input) integer scalar
</span><span class="comment">*</span><span class="comment">             These are the first and last rows and columns of an
</span><span class="comment">*</span><span class="comment">             isolated diagonal block upon which the QR sweep is to be
</span><span class="comment">*</span><span class="comment">             applied. It is assumed without a check that
</span><span class="comment">*</span><span class="comment">                       either KTOP = 1  or   H(KTOP,KTOP-1) = 0
</span><span class="comment">*</span><span class="comment">             and
</span><span class="comment">*</span><span class="comment">                       either KBOT = N  or   H(KBOT+1,KBOT) = 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      NSHFTS (input) integer scalar
</span><span class="comment">*</span><span class="comment">             NSHFTS gives the number of simultaneous shifts.  NSHFTS
</span><span class="comment">*</span><span class="comment">             must be positive and even.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      S      (input) COMPLEX array of size (NSHFTS)
</span><span class="comment">*</span><span class="comment">             S contains the shifts of origin that define the multi-
</span><span class="comment">*</span><span class="comment">             shift QR sweep.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      H      (input/output) COMPLEX array of size (LDH,N)
</span><span class="comment">*</span><span class="comment">             On input H contains a Hessenberg matrix.  On output a
</span><span class="comment">*</span><span class="comment">             multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied
</span><span class="comment">*</span><span class="comment">             to the isolated diagonal block in rows and columns KTOP
</span><span class="comment">*</span><span class="comment">             through KBOT.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      LDH    (input) integer scalar
</span><span class="comment">*</span><span class="comment">             LDH is the leading dimension of H just as declared in the
</span><span class="comment">*</span><span class="comment">             calling procedure.  LDH.GE.MAX(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      ILOZ   (input) INTEGER
</span><span class="comment">*</span><span class="comment">      IHIZ   (input) INTEGER
</span><span class="comment">*</span><span class="comment">             Specify the rows of Z to which transformations must be
</span><span class="comment">*</span><span class="comment">             applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      Z      (input/output) COMPLEX array of size (LDZ,IHI)
</span><span class="comment">*</span><span class="comment">             If WANTZ = .TRUE., then the QR Sweep unitary
</span><span class="comment">*</span><span class="comment">             similarity transformation is accumulated into
</span><span class="comment">*</span><span class="comment">             Z(ILOZ:IHIZ,ILO:IHI) from the right.
</span><span class="comment">*</span><span class="comment">             If WANTZ = .FALSE., then Z is unreferenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      LDZ    (input) integer scalar
</span><span class="comment">*</span><span class="comment">             LDA is the leading dimension of Z just as declared in
</span><span class="comment">*</span><span class="comment">             the calling procedure. LDZ.GE.N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      V      (workspace) COMPLEX array of size (LDV,NSHFTS/2)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      LDV    (input) integer scalar
</span><span class="comment">*</span><span class="comment">             LDV is the leading dimension of V as declared in the
</span><span class="comment">*</span><span class="comment">             calling procedure.  LDV.GE.3.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      U      (workspace) COMPLEX array of size
</span><span class="comment">*</span><span class="comment">             (LDU,3*NSHFTS-3)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      LDU    (input) integer scalar
</span><span class="comment">*</span><span class="comment">             LDU is the leading dimension of U just as declared in the
</span><span class="comment">*</span><span class="comment">             in the calling subroutine.  LDU.GE.3*NSHFTS-3.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      NH     (input) integer scalar
</span><span class="comment">*</span><span class="comment">             NH is the number of columns in array WH available for
</span><span class="comment">*</span><span class="comment">             workspace. NH.GE.1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      WH     (workspace) COMPLEX array of size (LDWH,NH)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      LDWH   (input) integer scalar
</span><span class="comment">*</span><span class="comment">             Leading dimension of WH just as declared in the
</span><span class="comment">*</span><span class="comment">             calling procedure.  LDWH.GE.3*NSHFTS-3.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      NV     (input) integer scalar
</span><span class="comment">*</span><span class="comment">             NV is the number of rows in WV agailable for workspace.
</span><span class="comment">*</span><span class="comment">             NV.GE.1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      WV     (workspace) COMPLEX array of size
</span><span class="comment">*</span><span class="comment">             (LDWV,3*NSHFTS-3)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      LDWV   (input) integer scalar
</span><span class="comment">*</span><span class="comment">             LDWV is the leading dimension of WV as declared in the
</span><span class="comment">*</span><span class="comment">             in the calling subroutine.  LDWV.GE.NV.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ================================================================
</span><span class="comment">*</span><span class="comment">     Based on contributions by
</span><span class="comment">*</span><span class="comment">        Karen Braman and Ralph Byers, Department of Mathematics,
</span><span class="comment">*</span><span class="comment">        University of Kansas, USA
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ============================================================
</span><span class="comment">*</span><span class="comment">     Reference:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
</span><span class="comment">*</span><span class="comment">     Algorithm Part I: Maintaining Well Focused Shifts, and
</span><span class="comment">*</span><span class="comment">     Level 3 Performance, SIAM Journal of Matrix Analysis,
</span><span class="comment">*</span><span class="comment">     volume 23, pages 929--947, 2002.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ============================================================
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      COMPLEX            ZERO, ONE
      PARAMETER          ( ZERO = ( 0.0e0, 0.0e0 ),
     $                   ONE = ( 1.0e0, 0.0e0 ) )
      REAL               RZERO, RONE
      PARAMETER          ( RZERO = 0.0e0, RONE = 1.0e0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      COMPLEX            ALPHA, BETA, CDUM, REFSUM
      REAL               H11, H12, H21, H22, SAFMAX, SAFMIN, SCL,
     $                   SMLNUM, TST1, TST2, ULP
      INTEGER            I2, I4, INCOL, J, J2, J4, JBOT, JCOL, JLEN,
     $                   JROW, JTOP, K, K1, KDU, KMS, KNZ, KRCOL, KZS,
     $                   M, M22, MBOT, MEND, MSTART, MTOP, NBMPS, NDCOL,
     $                   NS, NU
      LOGICAL            ACCUM, BLK22, BMP22
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      REAL               <a name="SLAMCH.156"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
      EXTERNAL           <a name="SLAMCH.157"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span><span class="comment">*</span><span class="comment">
</span>      INTRINSIC          ABS, AIMAG, CONJG, MAX, MIN, MOD, REAL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      COMPLEX            VT( 3 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           CGEMM, <a name="CLACPY.167"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>, <a name="CLAQR1.167"></a><a href="claqr1.f.html#CLAQR1.1">CLAQR1</a>, <a name="CLARFG.167"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</a>, <a name="CLASET.167"></a><a href="claset.f.html#CLASET.1">CLASET</a>, CTRMM,
     $                   <a name="SLABAD.168"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Functions ..
</span>      REAL               CABS1
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Function definitions ..
</span>      CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== If there are no shifts, then there is nothing to do. ====
</span><span class="comment">*</span><span class="comment">
</span>      IF( NSHFTS.LT.2 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== If the active block is empty or 1-by-1, then there
</span><span class="comment">*</span><span class="comment">     .    is nothing to do. ====
</span><span class="comment">*</span><span class="comment">
</span>      IF( KTOP.GE.KBOT )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== NSHFTS is supposed to be even, but if is odd,
</span><span class="comment">*</span><span class="comment">     .    then simply reduce it by one.  ====
</span><span class="comment">*</span><span class="comment">
</span>      NS = NSHFTS - MOD( NSHFTS, 2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Machine constants for deflation ====
</span><span class="comment">*</span><span class="comment">
</span>      SAFMIN = <a name="SLAMCH.196"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'SAFE MINIMUM'</span> )
      SAFMAX = RONE / SAFMIN
      CALL <a name="SLABAD.198"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>( SAFMIN, SAFMAX )
      ULP = <a name="SLAMCH.199"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'PRECISION'</span> )
      SMLNUM = SAFMIN*( REAL( N ) / ULP )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Use accumulated reflections to update far-from-diagonal
</span><span class="comment">*</span><span class="comment">     .    entries ? ====
</span><span class="comment">*</span><span class="comment">
</span>      ACCUM = ( KACC22.EQ.1 ) .OR. ( KACC22.EQ.2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== If so, exploit the 2-by-2 block structure? ====
</span><span class="comment">*</span><span class="comment">
</span>      BLK22 = ( NS.GT.2 ) .AND. ( KACC22.EQ.2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== clear trash ====
</span><span class="comment">*</span><span class="comment">
</span>      IF( KTOP+2.LE.KBOT )
     $   H( KTOP+2, KTOP ) = ZERO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== NBMPS = number of 2-shift bulges in the chain ====
</span><span class="comment">*</span><span class="comment">
</span>      NBMPS = NS / 2
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== KDU = width of slab ====
</span><span class="comment">*</span><span class="comment">
</span>      KDU = 6*NBMPS - 3
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== Create and chase chains of NBMPS bulges ====
</span><span class="comment">*</span><span class="comment">
</span>      DO 210 INCOL = 3*( 1-NBMPS ) + KTOP - 1, KBOT - 2, 3*NBMPS - 2
         NDCOL = INCOL + KDU
         IF( ACCUM )
     $      CALL <a name="CLASET.229"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'ALL'</span>, KDU, KDU, ZERO, ONE, U, LDU )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Near-the-diagonal bulge chase.  The following loop
</span><span class="comment">*</span><span class="comment">        .    performs the near-the-diagonal part of a small bulge
</span><span class="comment">*</span><span class="comment">        .    multi-shift QR sweep.  Each 6*NBMPS-2 column diagonal
</span><span class="comment">*</span><span class="comment">        .    chunk extends from column INCOL to column NDCOL
</span><span class="comment">*</span><span class="comment">        .    (including both column INCOL and column NDCOL). The
</span><span class="comment">*</span><span class="comment">        .    following loop chases a 3*NBMPS column long chain of
</span><span class="comment">*</span><span class="comment">        .    NBMPS bulges 3*NBMPS-2 columns to the right.  (INCOL
</span><span class="comment">*</span><span class="comment">        .    may be less than KTOP and and NDCOL may be greater than
</span><span class="comment">*</span><span class="comment">        .    KBOT indicating phantom columns from which to chase
</span><span class="comment">*</span><span class="comment">        .    bulges before they are actually introduced or to which
</span><span class="comment">*</span><span class="comment">        .    to chase bulges beyond column KBOT.)  ====
</span><span class="comment">*</span><span class="comment">
</span>         DO 140 KRCOL = INCOL, MIN( INCOL+3*NBMPS-3, KBOT-2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Bulges number MTOP to MBOT are active double implicit
</span><span class="comment">*</span><span class="comment">           .    shift bulges.  There may or may not also be small
</span><span class="comment">*</span><span class="comment">           .    2-by-2 bulge, if there is room.  The inactive bulges
</span><span class="comment">*</span><span class="comment">           .    (if any) must wait until the active bulges have moved
</span><span class="comment">*</span><span class="comment">           .    down the diagonal to make room.  The phantom matrix
</span><span class="comment">*</span><span class="comment">           .    paradigm described above helps keep track.  ====
</span><span class="comment">*</span><span class="comment">
</span>            MTOP = MAX( 1, ( ( KTOP-1 )-KRCOL+2 ) / 3+1 )
            MBOT = MIN( NBMPS, ( KBOT-KRCOL ) / 3 )
            M22 = MBOT + 1
            BMP22 = ( MBOT.LT.NBMPS ) .AND. ( KRCOL+3*( M22-1 ) ).EQ.
     $              ( KBOT-2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Generate reflections to chase the chain right
</span><span class="comment">*</span><span class="comment">           .    one column.  (The minimum value of K is KTOP-1.) ====
</span><span class="comment">*</span><span class="comment">
</span>            DO 10 M = MTOP, MBOT
               K = KRCOL + 3*( M-1 )
               IF( K.EQ.KTOP-1 ) THEN
                  CALL <a name="CLAQR1.264"></a><a href="claqr1.f.html#CLAQR1.1">CLAQR1</a>( 3, H( KTOP, KTOP ), LDH, S( 2*M-1 ),
     $                         S( 2*M ), V( 1, M ) )
                  ALPHA = V( 1, M )
                  CALL <a name="CLARFG.267"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</a>( 3, ALPHA, V( 2, M ), 1, V( 1, M ) )
               ELSE
                  BETA = H( K+1, K )
                  V( 2, M ) = H( K+2, K )
                  V( 3, M ) = H( K+3, K )
                  CALL <a name="CLARFG.272"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</a>( 3, BETA, V( 2, M ), 1, V( 1, M ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== A Bulge may collapse because of vigilant
</span><span class="comment">*</span><span class="comment">                 .    deflation or destructive underflow.  (The
</span><span class="comment">*</span><span class="comment">                 .    initial bulge is always collapsed.) Use
</span><span class="comment">*</span><span class="comment">                 .    the two-small-subdiagonals trick to try
</span><span class="comment">*</span><span class="comment">                 .    to get it started again. If V(2,M).NE.0 and
</span><span class="comment">*</span><span class="comment">                 .    V(3,M) = H(K+3,K+1) = H(K+3,K+2) = 0, then
</span><span class="comment">*</span><span class="comment">                 .    this bulge is collapsing into a zero
</span><span class="comment">*</span><span class="comment">                 .    subdiagonal.  It will be restarted next
</span><span class="comment">*</span><span class="comment">                 .    trip through the loop.)
</span><span class="comment">*</span><span class="comment">
</span>                  IF( V( 1, M ).NE.ZERO .AND.
     $                ( V( 3, M ).NE.ZERO .OR. ( H( K+3,
     $                K+1 ).EQ.ZERO .AND. H( K+3, K+2 ).EQ.ZERO ) ) )
     $                 THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    ==== Typical case: not collapsed (yet). ====
</span><span class="comment">*</span><span class="comment">
</span>                     H( K+1, K ) = BETA
                     H( K+2, K ) = ZERO
                     H( K+3, K ) = ZERO
                  ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    ==== Atypical case: collapsed.  Attempt to
</span><span class="comment">*</span><span class="comment">                    .    reintroduce ignoring H(K+1,K).  If the
</span><span class="comment">*</span><span class="comment">                    .    fill resulting from the new reflector
</span><span class="comment">*</span><span class="comment">                    .    is too large, then abandon it.
</span><span class="comment">*</span><span class="comment">                    .    Otherwise, use the new one. ====
</span><span class="comment">*</span><span class="comment">
</span>                     CALL <a name="CLAQR1.302"></a><a href="claqr1.f.html#CLAQR1.1">CLAQR1</a>( 3, H( K+1, K+1 ), LDH, S( 2*M-1 ),
     $                            S( 2*M ), VT )
                     SCL = CABS1( VT( 1 ) ) + CABS1( VT( 2 ) ) +
     $                     CABS1( VT( 3 ) )
                     IF( SCL.NE.RZERO ) THEN
                        VT( 1 ) = VT( 1 ) / SCL
                        VT( 2 ) = VT( 2 ) / SCL
                        VT( 3 ) = VT( 3 ) / SCL
                     END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    ==== The following is the traditional and
</span><span class="comment">*</span><span class="comment">                    .    conservative two-small-subdiagonals
</span><span class="comment">*</span><span class="comment">                    .    test.  ====
</span><span class="comment">*</span><span class="comment">                    .
</span>                     IF( CABS1( H( K+1, K ) )*
     $                   ( CABS1( VT( 2 ) )+CABS1( VT( 3 ) ) ).GT.ULP*
     $                   CABS1( VT( 1 ) )*( CABS1( H( K,
     $                   K ) )+CABS1( H( K+1, K+1 ) )+CABS1( H( K+2,
     $                   K+2 ) ) ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       ==== Starting a new bulge here would
</span><span class="comment">*</span><span class="comment">                       .    create non-negligible fill.   If
</span><span class="comment">*</span><span class="comment">                       .    the old reflector is diagonal (only
</span><span class="comment">*</span><span class="comment">                       .    possible with underflows), then
</span><span class="comment">*</span><span class="comment">                       .    change it to I.  Otherwise, use
</span><span class="comment">*</span><span class="comment">                       .    it with trepidation. ====
</span><span class="comment">*</span><span class="comment">
</span>                        IF( V( 2, M ).EQ.ZERO .AND. V( 3, M ).EQ.ZERO )
     $                       THEN
                           V( 1, M ) = ZERO
                        ELSE
                           H( K+1, K ) = BETA
                           H( K+2, K ) = ZERO
                           H( K+3, K ) = ZERO
                        END IF
                     ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       ==== Stating a new bulge here would
</span><span class="comment">*</span><span class="comment">                       .    create only negligible fill.
</span><span class="comment">*</span><span class="comment">                       .    Replace the old reflector with
</span><span class="comment">*</span><span class="comment">                       .    the new one. ====
</span><span class="comment">*</span><span class="comment">
</span>                        ALPHA = VT( 1 )
                        CALL <a name="CLARFG.345"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</a>( 3, ALPHA, VT( 2 ), 1, VT( 1 ) )
                        REFSUM = H( K+1, K ) +
     $                           H( K+2, K )*CONJG( VT( 2 ) ) +
     $                           H( K+3, K )*CONJG( VT( 3 ) )
                        H( K+1, K ) = H( K+1, K ) -
     $                                CONJG( VT( 1 ) )*REFSUM
                        H( K+2, K ) = ZERO
                        H( K+3, K ) = ZERO
                        V( 1, M ) = VT( 1 )
                        V( 2, M ) = VT( 2 )
                        V( 3, M ) = VT( 3 )
                     END IF
                  END IF
               END IF
   10       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Generate a 2-by-2 reflection, if needed. ====
</span><span class="comment">*</span><span class="comment">
</span>            K = KRCOL + 3*( M22-1 )
            IF( BMP22 ) THEN
               IF( K.EQ.KTOP-1 ) THEN
                  CALL <a name="CLAQR1.366"></a><a href="claqr1.f.html#CLAQR1.1">CLAQR1</a>( 2, H( K+1, K+1 ), LDH, S( 2*M22-1 ),
     $                         S( 2*M22 ), V( 1, M22 ) )
                  BETA = V( 1, M22 )
                  CALL <a name="CLARFG.369"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</a>( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) )
               ELSE
                  BETA = H( K+1, K )
                  V( 2, M22 ) = H( K+2, K )
                  CALL <a name="CLARFG.373"></a><a href="clarfg.f.html#CLARFG.1">CLARFG</a>( 2, BETA, V( 2, M22 ), 1, V( 1, M22 ) )
                  H( K+1, K ) = BETA
                  H( K+2, K ) = ZERO
               END IF
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Initialize V(1,M22) here to avoid possible undefined
</span><span class="comment">*</span><span class="comment">              .    variable problems later. ====
</span><span class="comment">*</span><span class="comment">
</span>               V( 1, M22 ) = ZERO
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Multiply H by reflections from the left ====
</span><span class="comment">*</span><span class="comment">
</span>            IF( ACCUM ) THEN
               JBOT = MIN( NDCOL, KBOT )
            ELSE IF( WANTT ) THEN
               JBOT = N
            ELSE
               JBOT = KBOT
            END IF
            DO 30 J = MAX( KTOP, KRCOL ), JBOT
               MEND = MIN( MBOT, ( J-KRCOL+2 ) / 3 )
               DO 20 M = MTOP, MEND
                  K = KRCOL + 3*( M-1 )
                  REFSUM = CONJG( V( 1, M ) )*
     $                     ( H( K+1, J )+CONJG( V( 2, M ) )*H( K+2, J )+
     $                     CONJG( V( 3, M ) )*H( K+3, J ) )
                  H( K+1, J ) = H( K+1, J ) - REFSUM
                  H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M )
                  H( K+3, J ) = H( K+3, J ) - REFSUM*V( 3, M )
   20          CONTINUE
   30       CONTINUE
            IF( BMP22 ) THEN
               K = KRCOL + 3*( M22-1 )
               DO 40 J = MAX( K+1, KTOP ), JBOT
                  REFSUM = CONJG( V( 1, M22 ) )*
     $                     ( H( K+1, J )+CONJG( V( 2, M22 ) )*
     $                     H( K+2, J ) )
                  H( K+1, J ) = H( K+1, J ) - REFSUM
                  H( K+2, J ) = H( K+2, J ) - REFSUM*V( 2, M22 )
   40          CONTINUE
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Multiply H by reflections from the right.
</span><span class="comment">*</span><span class="comment">           .    Delay filling in the last row until the
</span><span class="comment">*</span><span class="comment">           .    vigilant deflation check is complete. ====
</span><span class="comment">*</span><span class="comment">
</span>            IF( ACCUM ) THEN
               JTOP = MAX( KTOP, INCOL )
            ELSE IF( WANTT ) THEN
               JTOP = 1
            ELSE
               JTOP = KTOP
            END IF
            DO 80 M = MTOP, MBOT
               IF( V( 1, M ).NE.ZERO ) THEN
                  K = KRCOL + 3*( M-1 )
                  DO 50 J = JTOP, MIN( KBOT, K+3 )
                     REFSUM = V( 1, M )*( H( J, K+1 )+V( 2, M )*
     $                        H( J, K+2 )+V( 3, M )*H( J, K+3 ) )
                     H( J, K+1 ) = H( J, K+1 ) - REFSUM
                     H( J, K+2 ) = H( J, K+2 ) -
     $                             REFSUM*CONJG( V( 2, M ) )
                     H( J, K+3 ) = H( J, K+3 ) -
     $                             REFSUM*CONJG( V( 3, M ) )
   50             CONTINUE
<span class="comment">*</span><span class="comment">
</span>                  IF( ACCUM ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    ==== Accumulate U. (If necessary, update Z later
</span><span class="comment">*</span><span class="comment">                    .    with with an efficient matrix-matrix
</span><span class="comment">*</span><span class="comment">                    .    multiply.) ====
</span><span class="comment">*</span><span class="comment">
</span>                     KMS = K - INCOL
                     DO 60 J = MAX( 1, KTOP-INCOL ), KDU
                        REFSUM = V( 1, M )*( U( J, KMS+1 )+V( 2, M )*
     $                           U( J, KMS+2 )+V( 3, M )*U( J, KMS+3 ) )
                        U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM
                        U( J, KMS+2 ) = U( J, KMS+2 ) -
     $                                  REFSUM*CONJG( V( 2, M ) )
                        U( J, KMS+3 ) = U( J, KMS+3 ) -
     $                                  REFSUM*CONJG( V( 3, M ) )
   60                CONTINUE
                  ELSE IF( WANTZ ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    ==== U is not accumulated, so update Z
</span><span class="comment">*</span><span class="comment">                    .    now by multiplying by reflections
</span><span class="comment">*</span><span class="comment">                    .    from the right. ====
</span><span class="comment">*</span><span class="comment">
</span>                     DO 70 J = ILOZ, IHIZ
                        REFSUM = V( 1, M )*( Z( J, K+1 )+V( 2, M )*
     $                           Z( J, K+2 )+V( 3, M )*Z( J, K+3 ) )
                        Z( J, K+1 ) = Z( J, K+1 ) - REFSUM
                        Z( J, K+2 ) = Z( J, K+2 ) -
     $                                REFSUM*CONJG( V( 2, M ) )
                        Z( J, K+3 ) = Z( J, K+3 ) -
     $                                REFSUM*CONJG( V( 3, M ) )
   70                CONTINUE
                  END IF
               END IF
   80       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Special case: 2-by-2 reflection (if needed) ====
</span><span class="comment">*</span><span class="comment">
</span>            K = KRCOL + 3*( M22-1 )
            IF( BMP22 .AND. ( V( 1, M22 ).NE.ZERO ) ) THEN
               DO 90 J = JTOP, MIN( KBOT, K+3 )
                  REFSUM = V( 1, M22 )*( H( J, K+1 )+V( 2, M22 )*
     $                     H( J, K+2 ) )
                  H( J, K+1 ) = H( J, K+1 ) - REFSUM
                  H( J, K+2 ) = H( J, K+2 ) -
     $                          REFSUM*CONJG( V( 2, M22 ) )
   90          CONTINUE
<span class="comment">*</span><span class="comment">
</span>               IF( ACCUM ) THEN
                  KMS = K - INCOL
                  DO 100 J = MAX( 1, KTOP-INCOL ), KDU
                     REFSUM = V( 1, M22 )*( U( J, KMS+1 )+V( 2, M22 )*
     $                        U( J, KMS+2 ) )
                     U( J, KMS+1 ) = U( J, KMS+1 ) - REFSUM
                     U( J, KMS+2 ) = U( J, KMS+2 ) -
     $                               REFSUM*CONJG( V( 2, M22 ) )
  100             CONTINUE
               ELSE IF( WANTZ ) THEN
                  DO 110 J = ILOZ, IHIZ
                     REFSUM = V( 1, M22 )*( Z( J, K+1 )+V( 2, M22 )*
     $                        Z( J, K+2 ) )
                     Z( J, K+1 ) = Z( J, K+1 ) - REFSUM
                     Z( J, K+2 ) = Z( J, K+2 ) -
     $                             REFSUM*CONJG( V( 2, M22 ) )
  110             CONTINUE
               END IF
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Vigilant deflation check ====
</span><span class="comment">*</span><span class="comment">
</span>            MSTART = MTOP
            IF( KRCOL+3*( MSTART-1 ).LT.KTOP )
     $         MSTART = MSTART + 1
            MEND = MBOT
            IF( BMP22 )
     $         MEND = MEND + 1
            IF( KRCOL.EQ.KBOT-2 )
     $         MEND = MEND + 1
            DO 120 M = MSTART, MEND
               K = MIN( KBOT-1, KRCOL+3*( M-1 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== The following convergence test requires that
</span><span class="comment">*</span><span class="comment">              .    the tradition small-compared-to-nearby-diagonals
</span><span class="comment">*</span><span class="comment">              .    criterion and the Ahues &amp; Tisseur (LAWN 122, 1997)
</span><span class="comment">*</span><span class="comment">              .    criteria both be satisfied.  The latter improves
</span><span class="comment">*</span><span class="comment">              .    accuracy in some examples. Falling back on an
</span><span class="comment">*</span><span class="comment">              .    alternate convergence criterion when TST1 or TST2
</span><span class="comment">*</span><span class="comment">              .    is zero (as done here) is traditional but probably
</span><span class="comment">*</span><span class="comment">              .    unnecessary. ====
</span><span class="comment">*</span><span class="comment">
</span>               IF( H( K+1, K ).NE.ZERO ) THEN
                  TST1 = CABS1( H( K, K ) ) + CABS1( H( K+1, K+1 ) )
                  IF( TST1.EQ.RZERO ) THEN
                     IF( K.GE.KTOP+1 )
     $                  TST1 = TST1 + CABS1( H( K, K-1 ) )
                     IF( K.GE.KTOP+2 )
     $                  TST1 = TST1 + CABS1( H( K, K-2 ) )
                     IF( K.GE.KTOP+3 )
     $                  TST1 = TST1 + CABS1( H( K, K-3 ) )
                     IF( K.LE.KBOT-2 )
     $                  TST1 = TST1 + CABS1( H( K+2, K+1 ) )
                     IF( K.LE.KBOT-3 )
     $                  TST1 = TST1 + CABS1( H( K+3, K+1 ) )
                     IF( K.LE.KBOT-4 )
     $                  TST1 = TST1 + CABS1( H( K+4, K+1 ) )
                  END IF
                  IF( CABS1( H( K+1, K ) ).LE.MAX( SMLNUM, ULP*TST1 ) )
     $                 THEN
                     H12 = MAX( CABS1( H( K+1, K ) ),
     $                     CABS1( H( K, K+1 ) ) )
                     H21 = MIN( CABS1( H( K+1, K ) ),
     $                     CABS1( H( K, K+1 ) ) )
                     H11 = MAX( CABS1( H( K+1, K+1 ) ),
     $                     CABS1( H( K, K )-H( K+1, K+1 ) ) )
                     H22 = MIN( CABS1( H( K+1, K+1 ) ),
     $                     CABS1( H( K, K )-H( K+1, K+1 ) ) )
                     SCL = H11 + H12
                     TST2 = H22*( H11 / SCL )
<span class="comment">*</span><span class="comment">
</span>                     IF( TST2.EQ.RZERO .OR. H21*( H12 / SCL ).LE.
     $                   MAX( SMLNUM, ULP*TST2 ) )H( K+1, K ) = ZERO
                  END IF
               END IF
  120       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== Fill in the last row of each bulge. ====
</span><span class="comment">*</span><span class="comment">
</span>            MEND = MIN( NBMPS, ( KBOT-KRCOL-1 ) / 3 )
            DO 130 M = MTOP, MEND
               K = KRCOL + 3*( M-1 )
               REFSUM = V( 1, M )*V( 3, M )*H( K+4, K+3 )
               H( K+4, K+1 ) = -REFSUM
               H( K+4, K+2 ) = -REFSUM*CONJG( V( 2, M ) )
               H( K+4, K+3 ) = H( K+4, K+3 ) - REFSUM*CONJG( V( 3, M ) )
  130       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           ==== End of near-the-diagonal bulge chase. ====
</span><span class="comment">*</span><span class="comment">
</span>  140    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        ==== Use U (if accumulated) to update far-from-diagonal
</span><span class="comment">*</span><span class="comment">        .    entries in H.  If required, use U to update Z as
</span><span class="comment">*</span><span class="comment">        .    well. ====
</span><span class="comment">*</span><span class="comment">
</span>         IF( ACCUM ) THEN
            IF( WANTT ) THEN
               JTOP = 1
               JBOT = N
            ELSE
               JTOP = KTOP
               JBOT = KBOT
            END IF
            IF( ( .NOT.BLK22 ) .OR. ( INCOL.LT.KTOP ) .OR.
     $          ( NDCOL.GT.KBOT ) .OR. ( NS.LE.2 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Updates not exploiting the 2-by-2 block
</span><span class="comment">*</span><span class="comment">              .    structure of U.  K1 and NU keep track of
</span><span class="comment">*</span><span class="comment">              .    the location and size of U in the special
</span><span class="comment">*</span><span class="comment">              .    cases of introducing bulges and chasing
</span><span class="comment">*</span><span class="comment">              .    bulges off the bottom.  In these special
</span><span class="comment">*</span><span class="comment">              .    cases and in case the number of shifts
</span><span class="comment">*</span><span class="comment">              .    is NS = 2, there is no 2-by-2 block
</span><span class="comment">*</span><span class="comment">              .    structure to exploit.  ====
</span><span class="comment">*</span><span class="comment">
</span>               K1 = MAX( 1, KTOP-INCOL )
               NU = ( KDU-MAX( 0, NDCOL-KBOT ) ) - K1 + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Horizontal Multiply ====
</span><span class="comment">*</span><span class="comment">
</span>               DO 150 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH
                  JLEN = MIN( NH, JBOT-JCOL+1 )
                  CALL CGEMM( <span class="string">'C'</span>, <span class="string">'N'</span>, NU, JLEN, NU, ONE, U( K1, K1 ),
     $                        LDU, H( INCOL+K1, JCOL ), LDH, ZERO, WH,
     $                        LDWH )
                  CALL <a name="CLACPY.614"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'ALL'</span>, NU, JLEN, WH, LDWH,
     $                         H( INCOL+K1, JCOL ), LDH )
  150          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Vertical multiply ====
</span><span class="comment">*</span><span class="comment">
</span>               DO 160 JROW = JTOP, MAX( KTOP, INCOL ) - 1, NV
                  JLEN = MIN( NV, MAX( KTOP, INCOL )-JROW )
                  CALL CGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, JLEN, NU, NU, ONE,
     $                        H( JROW, INCOL+K1 ), LDH, U( K1, K1 ),
     $                        LDU, ZERO, WV, LDWV )
                  CALL <a name="CLACPY.625"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'ALL'</span>, JLEN, NU, WV, LDWV,
     $                         H( JROW, INCOL+K1 ), LDH )
  160          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Z multiply (also vertical) ====
</span><span class="comment">*</span><span class="comment">
</span>               IF( WANTZ ) THEN
                  DO 170 JROW = ILOZ, IHIZ, NV
                     JLEN = MIN( NV, IHIZ-JROW+1 )
                     CALL CGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, JLEN, NU, NU, ONE,
     $                           Z( JROW, INCOL+K1 ), LDZ, U( K1, K1 ),
     $                           LDU, ZERO, WV, LDWV )
                     CALL <a name="CLACPY.637"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'ALL'</span>, JLEN, NU, WV, LDWV,
     $                            Z( JROW, INCOL+K1 ), LDZ )
  170             CONTINUE
               END IF
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Updates exploiting U's 2-by-2 block structure.
</span><span class="comment">*</span><span class="comment">              .    (I2, I4, J2, J4 are the last rows and columns
</span><span class="comment">*</span><span class="comment">              .    of the blocks.) ====
</span><span class="comment">*</span><span class="comment">
</span>               I2 = ( KDU+1 ) / 2
               I4 = KDU
               J2 = I4 - I2
               J4 = KDU
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== KZS and KNZ deal with the band of zeros
</span><span class="comment">*</span><span class="comment">              .    along the diagonal of one of the triangular
</span><span class="comment">*</span><span class="comment">              .    blocks. ====
</span><span class="comment">*</span><span class="comment">
</span>               KZS = ( J4-J2 ) - ( NS+1 )
               KNZ = NS + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Horizontal multiply ====
</span><span class="comment">*</span><span class="comment">
</span>               DO 180 JCOL = MIN( NDCOL, KBOT ) + 1, JBOT, NH
                  JLEN = MIN( NH, JBOT-JCOL+1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== Copy bottom of H to top+KZS of scratch ====
</span><span class="comment">*</span><span class="comment">                  (The first KZS rows get multiplied by zero.) ====
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="CLACPY.667"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'ALL'</span>, KNZ, JLEN, H( INCOL+1+J2, JCOL ),
     $                         LDH, WH( KZS+1, 1 ), LDWH )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== Multiply by U21' ====
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="CLASET.672"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'ALL'</span>, KZS, JLEN, ZERO, ZERO, WH, LDWH )
                  CALL CTRMM( <span class="string">'L'</span>, <span class="string">'U'</span>, <span class="string">'C'</span>, <span class="string">'N'</span>, KNZ, JLEN, ONE,
     $                        U( J2+1, 1+KZS ), LDU, WH( KZS+1, 1 ),
     $                        LDWH )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== Multiply top of H by U11' ====
</span><span class="comment">*</span><span class="comment">
</span>                  CALL CGEMM( <span class="string">'C'</span>, <span class="string">'N'</span>, I2, JLEN, J2, ONE, U, LDU,
     $                        H( INCOL+1, JCOL ), LDH, ONE, WH, LDWH )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== Copy top of H bottom of WH ====
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="CLACPY.684"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'ALL'</span>, J2, JLEN, H( INCOL+1, JCOL ), LDH,
     $                         WH( I2+1, 1 ), LDWH )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== Multiply by U21' ====
</span><span class="comment">*</span><span class="comment">
</span>                  CALL CTRMM( <span class="string">'L'</span>, <span class="string">'L'</span>, <span class="string">'C'</span>, <span class="string">'N'</span>, J2, JLEN, ONE,
     $                        U( 1, I2+1 ), LDU, WH( I2+1, 1 ), LDWH )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== Multiply by U22 ====
</span><span class="comment">*</span><span class="comment">
</span>                  CALL CGEMM( <span class="string">'C'</span>, <span class="string">'N'</span>, I4-I2, JLEN, J4-J2, ONE,
     $                        U( J2+1, I2+1 ), LDU,
     $                        H( INCOL+1+J2, JCOL ), LDH, ONE,
     $                        WH( I2+1, 1 ), LDWH )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== Copy it back ====
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="CLACPY.701"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'ALL'</span>, KDU, JLEN, WH, LDWH,
     $                         H( INCOL+1, JCOL ), LDH )
  180          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Vertical multiply ====
</span><span class="comment">*</span><span class="comment">
</span>               DO 190 JROW = JTOP, MAX( INCOL, KTOP ) - 1, NV
                  JLEN = MIN( NV, MAX( INCOL, KTOP )-JROW )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== Copy right of H to scratch (the first KZS
</span><span class="comment">*</span><span class="comment">                 .    columns get multiplied by zero) ====
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="CLACPY.713"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'ALL'</span>, JLEN, KNZ, H( JROW, INCOL+1+J2 ),
     $                         LDH, WV( 1, 1+KZS ), LDWV )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== Multiply by U21 ====
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="CLASET.718"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'ALL'</span>, JLEN, KZS, ZERO, ZERO, WV, LDWV )
                  CALL CTRMM( <span class="string">'R'</span>, <span class="string">'U'</span>, <span class="string">'N'</span>, <span class="string">'N'</span>, JLEN, KNZ, ONE,
     $                        U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ),
     $                        LDWV )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== Multiply by U11 ====
</span><span class="comment">*</span><span class="comment">
</span>                  CALL CGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, JLEN, I2, J2, ONE,
     $                        H( JROW, INCOL+1 ), LDH, U, LDU, ONE, WV,
     $                        LDWV )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== Copy left of H to right of scratch ====
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="CLACPY.731"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'ALL'</span>, JLEN, J2, H( JROW, INCOL+1 ), LDH,
     $                         WV( 1, 1+I2 ), LDWV )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== Multiply by U21 ====
</span><span class="comment">*</span><span class="comment">
</span>                  CALL CTRMM( <span class="string">'R'</span>, <span class="string">'L'</span>, <span class="string">'N'</span>, <span class="string">'N'</span>, JLEN, I4-I2, ONE,
     $                        U( 1, I2+1 ), LDU, WV( 1, 1+I2 ), LDWV )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== Multiply by U22 ====
</span><span class="comment">*</span><span class="comment">
</span>                  CALL CGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, JLEN, I4-I2, J4-J2, ONE,
     $                        H( JROW, INCOL+1+J2 ), LDH,
     $                        U( J2+1, I2+1 ), LDU, ONE, WV( 1, 1+I2 ),
     $                        LDWV )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 ==== Copy it back ====
</span><span class="comment">*</span><span class="comment">
</span>                  CALL <a name="CLACPY.748"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'ALL'</span>, JLEN, KDU, WV, LDWV,
     $                         H( JROW, INCOL+1 ), LDH )
  190          CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              ==== Multiply Z (also vertical) ====
</span><span class="comment">*</span><span class="comment">
</span>               IF( WANTZ ) THEN
                  DO 200 JROW = ILOZ, IHIZ, NV
                     JLEN = MIN( NV, IHIZ-JROW+1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    ==== Copy right of Z to left of scratch (first
</span><span class="comment">*</span><span class="comment">                    .     KZS columns get multiplied by zero) ====
</span><span class="comment">*</span><span class="comment">
</span>                     CALL <a name="CLACPY.761"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'ALL'</span>, JLEN, KNZ,
     $                            Z( JROW, INCOL+1+J2 ), LDZ,
     $                            WV( 1, 1+KZS ), LDWV )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    ==== Multiply by U12 ====
</span><span class="comment">*</span><span class="comment">
</span>                     CALL <a name="CLASET.767"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'ALL'</span>, JLEN, KZS, ZERO, ZERO, WV,
     $                            LDWV )
                     CALL CTRMM( <span class="string">'R'</span>, <span class="string">'U'</span>, <span class="string">'N'</span>, <span class="string">'N'</span>, JLEN, KNZ, ONE,
     $                           U( J2+1, 1+KZS ), LDU, WV( 1, 1+KZS ),
     $                           LDWV )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    ==== Multiply by U11 ====
</span><span class="comment">*</span><span class="comment">
</span>                     CALL CGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, JLEN, I2, J2, ONE,
     $                           Z( JROW, INCOL+1 ), LDZ, U, LDU, ONE,
     $                           WV, LDWV )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    ==== Copy left of Z to right of scratch ====
</span><span class="comment">*</span><span class="comment">
</span>                     CALL <a name="CLACPY.781"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'ALL'</span>, JLEN, J2, Z( JROW, INCOL+1 ),
     $                            LDZ, WV( 1, 1+I2 ), LDWV )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    ==== Multiply by U21 ====
</span><span class="comment">*</span><span class="comment">
</span>                     CALL CTRMM( <span class="string">'R'</span>, <span class="string">'L'</span>, <span class="string">'N'</span>, <span class="string">'N'</span>, JLEN, I4-I2, ONE,
     $                           U( 1, I2+1 ), LDU, WV( 1, 1+I2 ),
     $                           LDWV )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    ==== Multiply by U22 ====
</span><span class="comment">*</span><span class="comment">
</span>                     CALL CGEMM( <span class="string">'N'</span>, <span class="string">'N'</span>, JLEN, I4-I2, J4-J2, ONE,
     $                           Z( JROW, INCOL+1+J2 ), LDZ,
     $                           U( J2+1, I2+1 ), LDU, ONE,
     $                           WV( 1, 1+I2 ), LDWV )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    ==== Copy the result back to Z ====
</span><span class="comment">*</span><span class="comment">
</span>                     CALL <a name="CLACPY.799"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'ALL'</span>, JLEN, KDU, WV, LDWV,
     $                            Z( JROW, INCOL+1 ), LDZ )
  200             CONTINUE
               END IF
            END IF
         END IF
  210 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ==== End of <a name="CLAQR5.807"></a><a href="claqr5.f.html#CLAQR5.1">CLAQR5</a> ====
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
